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ABSTRACT 

A new tuning-free subgrid-scale model, termed 'locally-averaged scale-dependent dynamic' (LASDD) 
model, is developed and implemented in large-eddy simulations (LESs) of stable boundary layers. The new 
model dynamically computes the Smagorinsky coefficient and the subgrid-scale Prandtl number based on the 
local dynamics of the resolved velocity and temperature fields. Overall, the agreement between the statistics 
of the LES-generated turbulence and some well-established empirical formulations and theoretical predictions 
(e.g., Nieuwstadt's local scaling hypothesis) is remarkable. The results show clear improvements over most of 
the traditional subgrid-scale models in the surface layer. Moreover, in contrast to previous large-eddy simulations 
of stable boundary layers that have strong dependence on grid resolution, the simulated statistics obtained with 
the LASDD model show relatively little resolution dependence for the range of grid sizes considered here. In 
essence, we show that the new LASDD model is a robust subgrid-scale parameterization for reliable, tuning-free 
simulations of stable boundary layers, even with relatively coarse resolutions. 



1. Introduction 

Atmospheric boundary layers (ABLs) are usually classi- 
fied into three types: neutral, convective and stable, based 
on atmospheric stability (buoyancy effects) and the domi- 
nant mechanism of turbulence generation (Stull 1988; Arya 
2001). The boundary layer becomes stably stratified when- 
ever the underlying surface is colder than the air. Under this 
atmospheric condition, turbulence is generated by shear and 
destroyed by negative buoyancy and viscosity. Because of 
this competition between shear and buoyancy effects, the 
strength of turbulence in the stable boundary layer (SBL) 
is much weaker in comparison to the neutral and convec- 
tive boundary layers. As a result, the stable boundary layer 
is also much shallower and characterized by smaller eddy 
motions. Stable boundary layer turbulence has not received 
much attention despite its scientifically intriguing nature and 
practical significance (e.g., numerical weather prediction - 
NWP, and pollutant transport). This might be attributed to 
the lack of adequate field or laboratory measurements, to 
the inevitable difficulties in numerical simulations (arising 
from small scales of motion due to stratification) and to the 
intrinsic complexity in its dynamics (e.g., occurrences of 
intermittency, Kelvin-Helmholtz instability, gravity waves, 
low-level jets, meandering motions etc.) (Hunt et al. 1996; 
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Mahrt 1998; Derbyshire 1999). Not surprisingly, today, 
there is a general consensus among researchers that our un- 
derstanding of the stable boundary layer (especially the very 
stable regime) is quite poor (Mahrt 1998; Derbyshire 1999; 
Holtslag 2003) and 'even small future advances justify more 
work' Pahrt 1998). 

In order to improve our understanding of SBL turbulence 
and to explore some of its inherent characteristics, in this 
study we make use of a contemporary numerical modeling 
approach, known as large-eddy simulation. Following the 
pioneering works of Deardorff Peardorff 1970, 1972, 1974, 1980) 
over the years LES has become an indispensable tool to 
study the ABL (e.g., Moeng 1984; Nieuwstadt et al. 1991; 
Andren et al. 1994; Mason 1994; Sullivan et al. 1994; 
Kosovic 1997; Albertson and Parlange 1999; Porte-Agel 
et al. 2000; Beare et al. 2004). However, until now LES 
models have not been sufficiently faithful in reproducing 
the characteristics of moderately and strongly stable atmo- 
spheric boundary layers (Saiki et al. 2000; Holtslag 2003). 
The main weakness of LES is associated with our limited 
ability to accurately account for the dynamics that are not 
explicitly resolved in the simulations. Under stable condi- 
tions - due to flow stratification - the characteristic size of 
the eddies becomes increasingly smaller with increasing at- 
mospheric stability, which eventually imposes an additional 
burden on the LES subgrid-scale (SGS) models. The re- 
cent GABLS (Global Energy and Water Cycle Experiment 
Atmospheric Boundary Layer Study) LES intercomparison 
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study (Beare et al. 2004) highlights that LESs of moderately 
stable BLs are quite sensitive to SGS models at a relatively 
fine resolution of 6.25 m. At a coarser resolution (12.5 m), 
a couple of commonly used SGS models even laminarized 
spuriously. Occasionally, laminarization was manifested by 
a near-linear (without any curvature) temperature profile; at 
times the SGS contributions to the total momentum or heat 
fluxes were larger than fifty percent in the interior of the 
boundary layer. This breakdown of traditional SGS mod- 
els undoubtedly calls for improved SGS parameterizations 
in order to make LES a more reliable tool to study stable 
boundary layers. The present study is devoted towards this 
goal. 

The organization of this paper is as follows. In Section 
2, we briefly describe the basic philosophy of large-eddy 
simulation. The newly developed locally-averaged scale- 
dependent dynamic (LASDD) SGS modeling approach is 
presented in Section 3. Simulations of stably stratified atmo- 
spheric boundary layers are presented in Section 4. Lastly, 
in Section 5, we summarize our research and elaborate on 
the prospects for future research in this subject area. 

2. Subgrid-Scale Modeling and SGS Parameter 
Estimation 

In rotation-influenced ABLs, the equations governing the 
conservation of momentum (using the Boussinesq approxi- 
mation) and temperature are: 
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where t is time, Xj is the spatial coordinate in the j-direction, 
uj is the velocity component in that direction, 9 is poten- 
tial temperature, 8 is reference surface potential temper- 
ature, p is dynamic pressure, 8^ is the Kronecker delta, 
^•3 is the alternating unit tensor, g is the gravitational ac- 
celeration, f c is the Coriolis parameter and Fi is a forcing 
term (e.g., Geostrophic wind or imposed mean pressure gra- 
dient). Molecular dissipation and diffusion have been ne- 
glected since the Reynolds number of the ABL is very high 
and no near-ground viscous processes are resolved. The ( ) 
is used to define a horizontal plane average. The (• • •) de- 
notes a spatial filtering operation, using a filter of character- 
istic width Af. These filtered equations are now amenable 
for numerical solution (LES) on a grid of mesh size A 9 , con- 
siderably larger than the smallest scale of turbulent motion 
(the so-called Kolmogorov scale). It is common practice to 
use the ratio of filter- width to grid-spacing, Af/A g = 1 or 
2 [see the Chapter 9 of Geurts (2003) for detailed discussion 
on this ratio and its impact on error dynamics]. In this study, 
we use a ratio of 2. 



The effects of the unresolved scales (smaller than A /) on 
the evolution of u.i and 9 appear in the SGS stress (see 
Equation 1) and the SGS flux qi (see Equation 2), respec- 
tively. They are defined as 



and 
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Note that the SGS stress and flux quantities are unknown and 
must be parameterized (using a SGS model) as a function of 
the resolved velocity and temperature fields. 

Eddy-viscosity models, the most popular SGS models, 
use the 'gradient hypothesis' and formulate the i j-component 
of the SGS stress tensor (the deviatoric part) as follows 
(Smagorinsky 1963; Geurts 2003): 
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where 
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is the resolved strain rate tensor and v t denotes the eddy- 
viscosity. It is well-known that the eddy-viscosity SGS mod- 
els give a poor prediction of the SGS stresses on a local 
level (see Sarghini et al. 1999; Meneveau and Katz 2000; 
Geurts 2003.) Their underlying assumption of strain rates 
being aligned with the SGS stress tensor is unrealistic (see 
Higgins et al. 2003 and the references therein). Further- 
more, these SGS models are purely dissipative, i.e., they do 
not allow local reverse energy transfer (known as 'backscat- 
ter') (Sarghini et al. 1999). Despite all these deficiencies, 
without any doubt, the eddy-viscosity models are the most 
commonly used SGS models in the atmospheric boundary 
layer community. 

From dimensional analysis, the eddy-viscosity (v t ) can be 
interpreted as the product of a characteristic velocity scale 
and a characteristic length scale (Geurts 2003). Different 
eddy-viscosity formulations basically use different velocity 
and length scales. The most popular eddy-viscosity formu- 
lation is the Smagorinsky model (Smagorinsky 1963): 



v t = (C s A f ) 2 



(6) 



where C5 is the so-called Smagorinsky coefficient, which 
is adjusted empirically or dynamically to account for shear, 
stratification and grid-resolution, and 



\S\ 



1/2 



is the magnitude of the resolved strain rate tensor. In con- 
trast to the Smagorinsky-type eddy-viscosity model, the tur- 
bulence kinetic energy (TKE) based eddy-viscosity model 
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utilizes (Moeng 1984; Sullivan et al. 1994; Sullivan et al. 
2003): 



(7) 



where Ck is a modeling coefficient, I is a length scale and 
Esgs is the SGS turbulence kinetic energy. This model- 
ing approach involves solving an extra prognostic equation 
for the SGS TKE. Based on the Kolmogorov's scaling laws, 
Wong and Lilly (1994) proposed yet another eddy-viscosity 
model: 

v t = C^Afe 1 / 3 = C e Ay\ (8) 

where e is the dissipation rate of energy and C £ is a model 
coefficient to be specified (or determined dynamically). There 
are numerous other formulations for eddy-viscosity existing 
in the literature (e.g., the Structure Function model of Metais 
and Lesieur 1992). An extensive review of these formula- 
tions is given by Sagaut (2001). 

Similar to the SGS stresses, the SGS heat fluxes are mod- 
eled with the eddy-diffusivity models as: 
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where Prscs is the SGS Prandtl number. 

The values of the Smagorinsky-type SGS model param- 
eters Cs and PrsGS are well established for homogeneous, 
isotropic turbulence (Lilly 1967). However, the value of 
Cs is expected to decrease with increasing mean shear and 
stratification. This has been confirmed by various recent 
field studies (Porte- Agel et al. 2001; Kleissl et al. 2003; 
Sullivan et al. 2003; Kleissl et al. 2004). In order to ac- 
count for this, application of the traditional eddy-viscosity 
model in LES of ABL flows (with strong shear near the 
ground and temperature-driven stratification) has tradition- 
ally involved the use of various types of wall-damping func- 
tions and stability corrections, which are either based on the 
phenomenological theory of turbulence or empirically de- 
rived from observational data (Mason 1994). For example, 
recently, based on the HATS (Horizontal Array Turbulence 
Study) field campaign data Kleissl et al. (2003, 2004) pro- 
posed the following empirical form for Cs- 
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where L is the Obukhov length, k is the von Karman con- 
stant, R is the ramp function, n = 3 and cq s=a 0.135. 
Another example would be Smagorinsky-type SGS models 
that impose both the wall-damping and stability corrections 
based on the Kansas field experiment data (see Brown et al. 
1994; Mason 1994; Mason and Brown 1999; Beare and 
MacVean 2004). In the case of TKE-based eddy-viscosity 
models, the length scale I is usually set equal to the filter 
width A f for unstable and neutral stratifications and equal to 



Ciy/Escs/N for stable stratification (Sullivan et al. 1994; 
Saiki et al. 2000; Sullivan et al. 2003), following the sugges- 
tion of Deardorff (1980). Here, C; is a coefficient to be pre- 
scribed and N is the Brunt-Vaisala frequency. When using 
this approach, the SGS model coefficients are often 'tuned' 
for different ABL flow conditions (Sullivan et al. 1994; Saiki 
et al. 2000; Sullivan et al. 2003). There also have been a 
few elegant attempts to derive shear and stability dependent 
length-scales directly from the phenomenological theory of 
turbulence (Hunt et al. 1988; Schumann 1991; Canuto and 
Cheng 1997). 'The adequacy of all these parameterizations 
for SGS fluxes remains relatively untested however' (Sulli- 
van et al. 2003). 

In the case of eddy-diffusivity SGS models, one needs 
to prescribe the stability dependence of the SGS Prandtl 
number (Prscs)- ln general, Prscs is found to increase 
under stable stratification, which is reflected in different 
SGS modeling approaches. For example, in a widely used 
Smagorinsky-type SGS model, the Prandtl number is in- 
creased from 0.44 in the free convection limit to 0.7 in neu- 
tral condition to 1 .0 in the very stable regime (Mason and 
Brown 1999). In contrast, the TKE-based SGS model uses 
(Moeng 1984; Sullivan et al. 1994; Saiki et al. 2000; Sul- 
livan et al. 2003): Pr SG s = A / /(A / + 21), where I is 
defined as before. This implies that the Pr S cs is 0.33 un- 
der convective and neutral conditions and varies from 0.33 
(weakly stable) to 1.0 (very stable) in the stably stratified 
regime. 

In summary, most of the conventional eddy-viscosity and 
eddy-diffusivity SGS modeling approaches involve param- 
eter tuning or a priori prescription in one way or another. 
An alternative approach is to use the 'dynamic' SGS model- 
ing approach of Germano (Germano et al. 1991; Germano 
1992; Lilly 1992). In this approach, one computes the value 
of the unknown SGS coefficient (e.g., the coefficient Cs in 
the Smagorinsky-type eddy-viscosity models) dynamically 
at every time and position in the flow. By looking at the 
dynamics of the flow at two different resolved scales and 
assuming scale similarity as well as scale invariance of the 
model coefficient, one can optimize its value. Thus, the dy- 
namic model avoids the need for a priori specification and 
tuning of the coefficient because it is evaluated directly from 
the resolved scales in an LES. In Moin et al. (1991), a similar 
dynamic procedure was applied to estimate the SGS scalar 
flux in compressible flows. In essence this procedure not 
only eliminates the need for any ad-hoc assumption about 
the SGS Prandtl number (Prscs) but also completely de- 
couples the SGS flux estimation from SGS stress computa- 
tion, which is highly desirable. 
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3. Locally-Averaged Scale-Dependent Dynamic 
Modeling Approach 

In the previous section, we mentioned that the SGS stress 
tensor (r^) at the filter scale (A/) is defined as: = 
HiUj — UiUj. In a seminal work, Germano (Germano et al. 
1991; Germano 1992) proposed to invoke an additional ex- 
plicit test filter of width a A f in order to dynamically com- 
pute the SGS coefficients. Consecutive filtering at scales Af 
and at a Af leads to a SGS turbulent stress tensor (Ty) at the 
test filter scale aA f. 



\ j — Ui Uj u 



where an overline (• • •) denotes filtering at a scale of aA/. 
From the definitions of r%j and Tij an algebraic relation can 
be formed, known in the literature as the Germano identity: 



~ Tij Tij. 



(12) 



This identity is then effectively used to dynamically ob- 
tain unknown SGS model coefficients. In the case of the 
Smagorinsky model, this identity yields 1 : 
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where My = 2A 2 .- 




(C 2 s) Af 



Mi. 



(13) 



Sij . If one 



assumes scale invariance, i.e., (C|) aA = {^s)a ^ er ~ 



manoetal. 1991), then the unknown coefficient (C|) A lmi 
be easily determined following the error minimization ap- 
proach of Lilly (1992): 



(<3)a, 
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Here the angular brackets (• • •) denote spatial averaging 
(Germano et al. 1991; Lilly 1992; Zang et al. 1993; Ghosal 
et al. 1995; Porte- Agel et al. 2000) or averaging over fluid 
pathlines (Meneveau et al. 1996). 

In a recent study, Porte- Agel et al. (2000) found that in a 
simulation of neutral boundary layer the scale-invariant dy- 
namic model is not dissipative enough in the near-ground 
region. It underestimates shear in that region and also yields 
excessively flat spectra (P orte-Agel et al. 2000). Moreover, 
the dynamically computed coefficient Cf show clear scale- 
dependence near the surface, i.e., ((7|) A ^ (^s) a A ■ 
Similar inferences are also obtained in the case of passive 
scalars (Porte-Agel 2004). Field observations by Kleissl 
et al. (2004) also support these results. 



1 Please note that here the variation of C| over the test filter scale has 
been implicitly neglected (Germano et al. 1991; Lilly 1992; Vreman et al. 
1994). 



This prompted Porte-Agel et al. (2000) to propose the 
scale-dependent dynamic SGS model. In this case, the scale- 

(c|) aA 

dependence parameter (3 — , T 1 is not assumed to be 

( C s) &f 

equal to one, rather it is determined dynamically. In order 
to implement this procedure, one needs to employ a second 
test filtering operation at a scale of a 2 Ay [denotedby (• • •)]. 
Invoking the Germano identity for the second time leads to: 
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This results in: 
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In the scale-dependent dynamic modeling approach, the fol- 
lowing assumption is made: 
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which is a much weaker assumption than the scale-invariance 
modeling assumption of [3 = 1. Now, from Equations 14 
and 16, one solves for the unknown parameter (3, which in 
turn is used to compute (C|) A/ utilizing Equation 14. Fur- 
ther technical details on this model could be found in Porte- 
Agel et al. (2000). In Porte-Agel (2004), similar formula- 
tions were derived for scalars. In this case, the lumped SGS 
coefficient CgPrgXa was determined dynamically and the 
scale-dependent parameter was termed as (3e- In the simu- 
lations of neutral boundary layers, the scale-dependent SGS 
model was found to exhibit 'right' dissipation behavior and 
more accurate spectra in the case of momentum Porte-Agel 
et al. 2000), as well as for passive scalars Porte-Agel 2004). 

In the present study we found that the original formula- 
tion of Porte-Agel et al. (2000) which involves (horizontal) 
planar averaging in Equation 14, suffers from an insufficient 
SGS dissipation problem in the outer layers in simulations 
of stable boundary layers. This could be attributed to decou- 
pling between horizontal planes under stratification. Inter- 
mittent, patchy turbulence in the strongly stable outer lay- 
ers might be another cause. This issue was resolved by us- 
ing a local formulation of the scale-dependent modeling ap- 
proach, named as the locally-averaged scale-dependent dy- 
namic (LASDD) model. The model coefficients (C| and 
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CgPrgQ S ) are obtained dynamically by averaging locally 
on the horizontal plane with a stencil of three by three grid 
points. Zang et al. (1993) followed a similar approach in 
the scale-invariant dynamic (i.e., f3 = 1) modeling of turbu- 
lent recirculating flows. To avoid numerical instabilities the 
coefficients C| and CgPr^^g are set to zero whenever the 
dynamic procedure yields negative values. This commonly 
used procedure is known as 'clipping' (Geurts 2003). The 
scale-dependence parameters {(3 and (3g) are determined dy- 
namically over horizontal planes to avoid the computational 
burden of computing them at every grid point in the flow. 
Solving for j3 or (3g involves a fifth-order polynomial. In- 
stead of the Newton-Raphson method used by Porte-Agel 
et al. (2000) and Porte-Agel (2004), we use a more robust 
eigenvalue based method (Press et al. 1992) to obtain the 
roots of this polynomial. In the infrequent events that an ap- 
propriate real root in the range of to 1 .2 is not found, we 
chose to invoke the scale-invariance assumption of (3 (or, f}g) 
= 1. 

Please note that, mathematically more rigorous (and com- 
putationally more expensive) local models are also proposed 
in the literature (Ghosal et al. 1995; Piomelli and Liu 1995; 
Meneveau et al. 1996). Their capabilities in the stably strat- 
ified atmospheric boundary layer simulations have yet to be 
tested. 
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erately stable boundary layers it led to development of un- 
physical profiles of various turbulent quantities (Saiki et al. 
2000). The failure was due to the collapse of SGS vertical 
heat flux near the surface. This prompted them to propose a 
two-part SGS model to represent the SGS heat flux similar 
to the SGS momentum flux model of Sullivan et al. (1994). 
However, even after the modifications, the simulations were 
found to be too sensitive to rapid cooling (Saiki et al. 2000). 

Kosovic and Curry (2000) simulated a clear-air, moder- 
ately stable boundary layer as it approaches quasi-steady 
state using the nonlinear SGS model of Kosovic (1997). 
Initial conditions consistent with the Beaufort Sea Arctic 
Stratus Experiment (BASE) observations were used. The 
first intercomparison study of LES-SGS models for stable 
boundary layer (Holtslag 2003; Beare et al. 2004), as part 
of the GABLS initiative, also used this case (slightly modi- 
fied). Eleven different models with very different SGS mod- 
eling options and different grid-resolutions (from 1 m to 
12.5 m) were run (Beare et al. 2004). In this paper, we 
also simulate this particular case with our newly proposed 
locally-averaged scale-dependent dynamic (LASDD) SGS 
model. We compare our results with theoretical predictions 
by Nieuwstadt (Nieuwstadt 1984a, b, 1985), field obser- 
vations, as well as with well-established empirical formu- 
lations. 



4. LES of Stably Stratified Boundary Layers 

The first LES of a stable boundary layer was performed 
by Mason and Derbyshire (1990). They used the traditional 
Smagorinsky-type SGS model with a constant SGS Prandtl 
number (Prscs — 0.5). Their results broadly supported the 
local scaling hypothesis of Nieuwstadt (Nieuwstadt 1984a; 
Nieuwstadt 1984b; Nieuwstadt 1985; Derbyshire 1990). 
However, one of their simulations showed the run-away 
cooling problem. The simulated surface temperature fell 
more than 30 K over 90 minutes. Brown et al. (1994) re- 
peated these simulations with their stochastic backscatter 
SGS model (with stability corrections). Their results defi- 
nitely showed some improvements (especially in the surface 
layer properties) when compared with the Mason and Der- 
byshire's simulations. 

Andren (1995) simulated weakly stable boundary layers 
with the TKE-based SGS model of Moeng (1984) and also 
with the two-part eddy-viscosity model developed by Sulli- 
van et al. (1994). The two-part eddy-viscosity model, which 
is a modified version of the TKE-based SGS model, was in 
better agreement with the surface-layer similarity theory. 

Recently, Saiki et al. (2000) attempted to simulate a mod- 
erately stable boundary layer with Sullivan et al.'s two-part 
eddy-viscosity model (Sullivan et al. 1994). Although, this 
particular SGS scheme has been previously used by Andren 
for weakly stable BLs (Andren 1995), in the case of mod- 



a. Description of the Simulations 

The GABLS LES intercomparison case study is described 
in detail in Beare et al. (2004). The boundary layer is 
driven by an imposed, uniform geostrophic wind (G = 8 
m s _1 ), with a surface cooling rate of 0.25 K per hour and 
attains a quasi-steady state in ~ 8-9 hours with a boundary 
layer depth of ~ 200 m. The initial mean potential tem- 
perature is 265 K upto 100 m with an overlying inversion 
of strength 0.01 K m _1 . The Coriolis parameter is set to 
f c = 1.39 x 10~ 4 s" 1 , corresponding to latitude 73° N. Our 
domain size is: (L x = L y = L z = 400 m). This domain 
is divided into: (1) N x x N y x N z = 32 x 32 x 32 nodes 
(i.e., = A y = A z = 12.5 m); (2) N x x N y x N z = 
64 x 64 x 64 nodes (i.e., A x = A y = A z = 6.25 
m); and (3) N x x N y x N z = 80 x 80 x 80 nodes (i.e., 
A x = A y = A z = 5 m). One of the objectives behind these 
simulations was to investigate the sensitivity of our results 
to grid resolution. 

A Galilean transformation is used to weaken the sta- 
bility constraints on the time step. The grid moves with 
(UGai, Voai) — (5.5, 0) m/s in the 32 3 and 64 3 nodes cases. 
In the case of 80 3 simulation, we have used (Uoai, Va a i) = 
(5, 0) m/s. The time steps (At) for our 32 3 , 64 3 , and 80 3 
simulations are 0.4, 0.2 and 0.14 seconds, respectively. 
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b. Description of the LES Code 

In this work, we have used a modified version of the LES 
code described in Albertson and Parlange (1999), Porte- 
Agel et al. (2000), and Porte-Agel (2004). The salient fea- 
tures of this code are as follows: 

• It solves the filtered Navier-Stokes equations written 
in rotational form (Orszag and Pao 1974). 

• Derivatives in the horizontal directions are computed 
using the Fourier Collocation method, while vertical 
derivatives are approximated with second-order cen- 
tral differences Canuto et al. 1988). 

• Dealiasing of the nonlinear terms in Fourier space is 
done using the 3/2 rule Canuto et al. 1988). 

• Explicit second-order Adams-Bashforth time advance- 
ment scheme is used C anut o et al. 1988). 

• Explicit spectral cutoff filtering is used. The ratio be- 
tween the filter-width (A/) and grid-spacing (A g ) is 
set to 2. Here a is taken to be equal to \f2. 

• Only Coriolis terms involving horizontal wind are 
considered. 

• Forcing is imposed by Geostrophic wind. 

• Staggered vertical grid is used. 

The lower boundary condition is based on the Monin- 
Obukhov similarity theory with a surface roughness length 
z . The instantaneous components of surface shear stresses 
t xz and T yz are represented as functions of the resolved ve- 
locities u and v, at the grid point immediately above the sur- 
face (i.e., at a height of z = A z /2 in our case): 



U{z) 
v{z)~ 



U(z) 



(18) 



(19) 



where is the friction velocity, which is computed from the 
mean horizontal wind speed U(z) = ((u 2 + v 2 ) 1 ^ 2 ) at the 
first vertical model level (z — A z /2) as follows: 

U(z)k 



h 2-' 



logf— 

In a similar manner, the surface heat flux is computed as: 

W*K [0 S — O(z)] 



(20) 



(wO s ) = 



log(£) + &„£ 



(21) 



TABLE 1 . Basic characteristics of the simulated SBLs during the 
last hour of simulation. 



Grid Points 


H{m) 


L(m) 


u r (ms i ) 


0, (K) 


32 x 32 x 32 


205 


113 


0.283 


0.047 


64 x 64 x 64 


185 


114 


0.276 


0.045 


80 x 80 x 80 


192 


122 


0.285 


0.045 



where 6 S and Q(z) denote the surface temperature and the 
mean resolved potential temperature at the first model level, 



respectively. Following the recommendations of the GABLS 
intercomparison study, the constants b m and bh were set to 
4.8 and 7.8, respectively. 

The upper boundary consists of a zero stress condition, 
whereas the lateral boundary condition assumes periodicity. 
A Rayleigh damping layer at 300 m is used following the 
GABLS case description. 

c. Results 

In this section we report the results of our tuning-free sim- 
ulations and attempt to evaluate them against the theoretical 
predictions and well established observations-based formu- 
lations. Our results show clear improvements over most of 
the traditional models in the surface layer. We would like 
to emphasize that in the surface layer the relative contribu- 
tion of the SGS to the overall turbulent fluxes is very large. 
This is also the location where gradients are much stronger. 
Hence, simulation results become very sensitive to SGS for- 
mulations near the ground. The comparisons we perform 
are extensive and address: (1) temporal evolution of simu- 
lated statistics, (2) various first order statistics of turbulent 
velocity and temperature fields, (3) second order statistics of 
turbulent velocity and temperature fields, and (4) character- 
istics of dynamically estimated SGS coefficients. Details of 
these comparisons are provided below. 

The boundary layer height (H), Obukhov length (L) and 
other characteristics of the simulated SBLs using the locally- 
averaged scale-dependent dynamic SGS model (averaged 
over the final hour of simulation) are given in Table 1 . Fol- 
lowing Kosovic and Curry (2000) and Beare et al. (2004), 
the boundary layer height (H) is defined as (1/0.95) times 
the height where the mean local stress falls to five percent of 
its surface value. From this table, it is apparent that the sim- 
ulated (bulk) boundary layer parameters are quite insensitive 
to the grid resolution. In LES this behavior is always desir- 
able and its existence is usually attributed to the strength of 
a SGS model. 

1) Temporal Evolution 

In Figures 1 and 2 the time series of surface momentum 
flux, surface buoyancy flux, Obukhov length (L) and bound- 
ary layer height (H) are shown. The surface momentum 
flux reaches quasi-equilibrium after 4 hours of simulation. 
On the other hand the Obukhov length and boundary layer 
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FIG. 1 . Time series of surface momentum flux (top) and surface 
buoyancy flux (bottom). 
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FIG. 2. Time series of Obukhov length (top) and boundary layer 
height (bottom). 



height equilibrate well before 2 hours of simulation. Since 
the surface boundary condition is prescribed by a fixed cool- 
ing rate rather than a fixed flux, it is anticipated that the sur- 
face buoyancy flux will keep on evolving with time. 

2) First-order Statistics 

The mean profiles of wind speed, wind angle and po- 
tential temperature averaged over the final hour (8-9 hours) 
of simulation, are shown in Figures 3 and 4. The super- 
geostrophic nocturnal jet near the top of the boundary layer, 
is in accordance with Nieuwstadt's theoretical model for 
'stationary' stable boundary layers [see Equation 17 of Nieuw- 
stadt (1985)]. However, the angle between the surface wind 
direction and the geostrophic wind simulated by our LES is 
~ 35 degrees. This is much smaller than Nieuwstadt's pre- 
diction of 60 degrees. The second-order closure model of 
Brost and Wyngaard (1978) also predicts a value of ~ 58 de- 



grees. In contrast, the results from the GABLS study Beare 
et al. 2004) and Kosovic and Curry (2000) are in agreement 
with our results. 

Nieuwstadt also derived the following mean temperature 
profile [Equation 21 of Nieuwstadt (1985)]: 



e-. 



0* 



Rig H 
uRi 2 L 



In 



('-*)■ 



(22) 



where Ri f and Ri g denote the flux and the gradient Richard- 
son numbers, respectively. 6* ^= — signifies the sur- 
face layer temperature scale. Equation 22 implies that the 
temperature profile exhibits positive curvature (<~ d 2 Q/dz 2 ), 
which is clearly visible in Figure 4. 

We would like to point out that Nieuwstadt's analytical 
model is based on the hypothesis that the gradient Richard- 
son number (Ri g ) and the flux Richardson number (Rif) 
are constant with height inside the stable boundary layer. 
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FIG. 3. Mean wind speed (top) and wind angle (bottom) profiles. 
These profiles are averaged over the last one hour of simulation. 



Nieuwstadt was aware of the fact that this hypothesis does 
not hold for the lower part of the boundary layer (Nieuw- 
stadt 1985). In fact, Ri g and Ri / should go to zero near the 
surface (Nieuwstadt 1985), as can be seen from our sim- 
ulations (Figure 5). The violation of the basic assumption 
in the proximity of the land surface might explain some of 
the discrepancies between the LES results and Nieuwstadt's 
predictions. 

The Richardson numbers represent the ratio of the amount 
of turbulent kinetic energy (TKE) destroyed by buoyancy 
forces to the amount of TKE generated by wind shear (Stull 
1988). The values of Rif are consistently higher than the 
corresponding Ri g values, which is expected (see below). 
In the interior part of the boundary layer, Ri g is more or less 
constant (~ 0.2), in accord with Nieuwstadt's assumption. 
However, Ri f increases monotonically and is higher than 
0.2 in the upper part of the boundary layer. The magnitudes 
of both these Richardson numbers increase sharply near the 
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FIG. 4. Mean temperature profiles. These profiles are averaged 
over the last one hour of simulation. 

top of the boundary layer and become more than 1 in the 
inversion layer. 

It is straightforward to show that the ratio between Ri g 
and Ri / is the turbulent Prandtl number (Pr t ) Perbyshire 
1999; Howell and Sun 1999): 

Km Ri Q 



Pr t = 



K 



H 



Ri, 



(23) 



where Km and Kh represent eddy diffusivities for momen- 
tum and heat flux, respectively. The dependence of Pr t on 
atmospheric stability is not strong Perbyshire 1999; How- 
ell and Sun 1999). Inside the boundary layer (up to <~ 150 
m), (almost) all our simulated results yield = Pr t ~ 0.7 
(not shown here). Based on phenomenological theories of 
turbulence Townsend (1976) and Yakhot and Orszag (1986) 
also derived Pr t = 0.7. However, in the surface layer our 
results show that the values of Pr t increase to ~ 1. This 
is consistent with 'Microfronts' field experimental data an- 
alyzed by Howell and Sun (1999). They found on average, 
the estimates of Prt at 3 m level are higher than at the 10 m 
level, indicating that the relative efficiency of turbulent mo- 
mentum transfer with respect to heat transfer increases in the 
proximity of the land surface Howell and Sun 1999). 

In SBL simulations, one can test the performance of a 
SGS model by plotting a local nondimensional shear: 



ML 




— \ ^~ + h- (24) 



(25) 



and nondimensional temperature gradient: 

kz (99 

as a function of local stability parameter (z/A) and com- 
paring with field-observations-based formulations. Here, A 
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Fig. 5. Mean profiles of gradient Richardson number (top) and 
flux Richardson number (bottom). These profiles are averaged over 
the last one hour of simulation. 



denotes the local Obukhov length. In this work, a subscript 
'l on the turbulence quantities (e.g., u*l) will be used to 
specify evalutation using local turbulence quantities - other- 
wise, surface values are implied. Recently, Mahrt and Vick- 
ers (2003) called this type of similarity theory 'hybrid sim- 
ilarity theory', since it approaches Monin-Obukhov similar- 
ity as z decreases and also conforms to z-less stratification 
as z — > oo. In Figure 6 we plot the 'hybrid' nondimension- 
alized gradients and compare them with the formulations by 
Businger et al. (1971): 



ML 



1 + 4.7^ 
A 



$i/L = 0.74 + 4.7-, 
and by Beljaars and Holtslag (1991): 
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A 



a + be 



(26) 
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(28) 



V 32 x 32 x 32 

n 64 x 64 x 64 

• 80 x 80 x 80 

Businger et al. 

- - Beljaars and Holtslag 




10 c 



V 32 x 32 x 32 

n 64 x 64 x 64 

• 80 x 80 x 80 

Businger et al. 

- - Beljaars and Holtslag 




z/A 



FIG. 6. Locally computed nondimensional gradients of velocity 
(top) and temperature (bottom) against the local stability parameter. 
These statistics are computed during the last one hour of simula- 
tion. The field-observations-based formulations given by Businger 
et al. (1971) and Beljaars and Holtslag (1991) are also shown. 



<S>HL 



,, 2az\ 1/2 , 
«(1+3 T 1 + * 



-d4 



1 + c-d 



A) 



(29) 

where the suggested values of the coefficients are (Bel- 
jaars and Holtslag 1991): a = 1,6 = 2/3, c = 5 and 
d = 0.35. Interestingly, both these simulated gradients plot- 
ted against z/A show slopes slightly smaller than the widely 
used Businger et al.'s formulations. Based on CASES99 
field observations data, Mahrt and Vickers (2003) found a 
slope of 3.7 [in contrast to 4.7 as proposed by Businger et al. 
(1971)], which also fits our LES results remarkably well. 
Previous studies, such as Beljaars and Holtslag (1991) also 
found that $ml and <&hl increase slower than Businger 
et al.'s formulations and they proposed the aforementioned 
nonlinear formulation (Beljaars and Holtslag 1991). 
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FIG. 7. Locally computed nondimensional gradients of velocity 
(top) and temperature (bottom) against local gradient Richardson 
number. These statistics are computed during the last one hour of 
simulation. The formulations given by Businger et al. (1971) are 
also shown. 



FIG. 8. Mean momentum flux (top) and buoyancy flux (bot- 
tom) profiles corresponding to the 80 3 simulation using the locally- 
averaged scale-dependent dynamic (LASDD) SGS model. These 
profiles are averaged over the last one hour of simulation. 



There exists another representation for the nondimen- 
sional gradients in terms of local gradient Richardson num- 
ber (Rig). Figure 7 once again shows that the agreement 
between our LES results and established formulations are 
quite satisfactory. In the literature, usually the critical gra- 
dient Richardson number (Ri gc ) is considered to be around 
0.25. For Ri gc > 0.25 turbulence is very weak (Stull 1988; 
Brown et al. 1994), which is also noticeable in Figure 7. 

3) Second-order Statistics 

The mean profiles of vertical momentum flux and buoy- 
ancy flux averaged over the final hour (8-9 hours) of sim- 
ulation are given in Figure 8. The dashed lines show the 
resolved flux, and the dotted lines denote the SGS contribu- 
tion to the flux. As would be anticipated, near the ground 



the SGS contribution is much larger than its resolved coun- 
terpart. In the GABLS intercomparison study, there is a sig- 
nificant spread between the total momentum and buoyancy 
flux profiles simulated with different models. In particular at 
the surface the mean momentum and buoyancy fluxes range 
from 0.06 to 0.08 m 2 s- 2 and -3.5 to -5.5 x 1CT 4 nrV 3 , re- 
spectively (Beare et al. 2004). Our simulated results also 
fall in these ranges. 

Perhaps more interesting is to explore the normalized flux 
profiles shown in Figure 9. Nieuwstadt's analytical model 
predictions are as follows (Nieuwstadt 1985): 

2 

^ = (1 - z/H f /2 (30) 

ui 

{ ^ = (l-z/H). (31) 
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Fig. 9. Mean normalized momentum flux (top) and normalized 
buoyancy flux (bottom) profiles from the simulations performed in 
the present work. These profiles are averaged over the last one 
hour of simulation. Theoretical predictions by Nieuwstadt (1985) 
are also shown for comparison. 



Our model simulated results are in close agreement with 
Nieuwstadt's predictions. 

Figures 10 and 11 show the (resolved) variances of ve- 
locity components and temperature. In the surface layer, 
the normalized resolved velocity variances are smaller than 
Nieuwstadt's field observations (Nieuwstadt 1984a, b). This 
is expected as the SGS contributions to these variances have 
not been added here. This will be done while studying 
Nieuwstadt's local scaling hypothesis and will be shown to 
be in excellent agreement with Nieuwstadt's observations 
(see below). We would like to point out that in contrast to our 
simulated results, the nonlinear SGS model simulations by 
Kosovic and Curry (2000) surprisingly yielded surface layer 
velocity variances (resolved plus SGS) which were <~ 40 
percent smaller than Nieuwstadt's observations. 
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Fig. 10. Resolved velocity variances from the simulations per- 
formed in the present work. These profiles are averaged over the 
last one hour of simulation. 



In his local scaling hypothesis, Nieuwstadt (Nieuwstadt 
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1984a, b, 1985) conjectured that under stable stratification, 
the local Obukhov length (A) based on local turbulent fluxes 
should be considered as a more fundamental length scale. 
Then, according to this hypothesis, dimensionless combina- 
tions of turbulent variables [gradients, fluxes, (co-)variances 
etc.] which are measured at the same height (z) could be 
expressed as 'universal' functions of the stability parameter, 
C (= z/A). Exact forms of these functions could be pre- 
dicted by dimensional analysis only in the asymptotic very 
stable case (( — > oo), as discussed below. In the very sta- 
ble regime (z-less condition), since any explicit dependence 
on z disappears, local scaling predicts that dimensionless 
turbulent quantities asymptotically approach constant values 
(Nieuwstadt 1984a, b, 1985). Local scaling could be viewed 
as a generalization of the well established Monin-Obukhov 
(M-O) similarity theory (Monin and Yaglom 1971; Sorbjan 
1989). M-O similarity theory is strictly valid in the surface 
layer (lowest 10% of the ABL), whereas local scaling de- 
scribes the turbulent structure of the entire SBL (Nieuwstadt 
1984a, b, 1985). 

Whether or not our LES-generated statistics support the 
local scaling hypothesis was studied extensively in our re- 
cent work (Basu et al. 2005). In that study, we also per- 
formed rigorous statistical analyses of field observations and 
wind-tunnel measurements in order to verify the validity of 
local scaling hypothesis under very stable conditions. An ex- 
tensive set of turbulence statistics, computed from field and 
wind-tunnel measurements and also from LES-generated 
datasets, supported the validity of the local scaling hypothe- 
sis (in the cases of traditional bottom-up as well as upside- 
down stable boundary layers over homogeneous, flat ter- 
rains). We demonstrated that non-turbulent effects need to 
be removed from field data while studying similarity hy- 



TABLE 2. Number of samples in each stability class. 
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32 a 
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S2 
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1 


2 


3 
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S4 


0.50-1.00 
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S5 


> 1.00 
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11 


14 



potheses, otherwise the results could be misleading (Basu 
et al. 2005). For completeness of the present paper, we de- 
cided to include some key local scaling results. For ease 
in representation, we categorize our LES-generated database 
based on local stabilities (z/A) (see Table 2). The class SI 
represents near neutral stability; while S5 corresponds to the 
very stable regime. 

In Figures 12 and 13 we plot the normalized standard 
deviation of velocity components and temperature, respec- 
tively. The results are presented using standard boxplot no- 
tation with marks at 95, 75, 50, 25, and 5 percentile of an 
empirical distribution. The SGS contributions to the total 
standard deviations are estimated following the approach of 
Mason (Mason 1989; Mason and Derbyshire 1990). 

It is quite evident from Figures 12 and 13 that the normal- 
ized standard deviation of the turbulence variables closely 
follows the local scaling predictions and also z-less strati- 
fication. In Table 3 we further report the median values of 
the turbulence statistics corresponding to the category S5. 
Loosely, these median values could be considered as the 
asymptotic z-less values, which are found to be remarkably 
close to Nieuwstadt's analytical predictions and also field 
observations (see Table 3). For an example, Nieuwstadt's 
theory predicts that the normalized vertical velocity stan- 
dard deviation asymptotically approaches <~ 1.4 in the z- 
less regime. In Basu et al. (2005), we observed this value to 
be in the narrow range of 1.4 to 1.6. Recently, Heinemann 
(2004) compiled a list (see Table 2 of their paper) of tur- 
bulence statistics under very stable conditions (( max <~ 25) 
reported by various researchers. They found an asymptotic 
value of ~ 1.6 for <t w /u*l, in accord with Sorbjan (1986). 

In Figure 14, we report the mutual correlations between 
u, w and 8. The z-less values are also reported in Table 
3. Once again, these values are very similar to the ones 
compiled by Heinemann (2004), our previous study (Basu 
et al. 2005), results of Sorbjan (1986) and theoretical predic- 
tions of Nieuwstadt (1984b). As a note, Kaimal and Finni- 
gan (1994) also report that for < ( < 1, r u g = 0.6, which 
is close to the values found in the present study. 

In light of the foregoing analyses it is certain that the lo- 
cal scaling hypothesis of Nieuwstadt, which has survived the 
last two decades, still holds for a wide range of stabilities and 
is well reproduced by our LES model. 
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FIG. 12. Normalized longitudinal (top), transverse (middle), and 
vertical (bottom) velocity standard deviations. These statistics are 
averaged over the last one hour of simulation. 

4) SGS Coefficients 

Figures 15 and 16 show the SGS coefficients: Cs, PrsGS, 
and the averaged scale-dependent parameters: 0, fig, dynam- 
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Fig. 13. Normalized temperature standard deviations. This 
statistic is averaged over the last one hour of simulation. 



ically obtained using the locally-averaged scale-dependent 
dynamic model. and 0$ are found to be significantly 
smaller than 1 in the entire boundary layer. This stresses the 
fact that the assumption of scale-invariance in anisotropic 
stably stratified flows is not appropriate. Indeed, both Cs 
and Pr S GS are found to be scale dependent. The scale- 
dependent parameters are also expected to depend on local 
stability as they decrease significantly in the strongly strati- 
fied inversion layer. Cs is found to decrease with increasing 
atmospheric stability, consistent with recent field observa- 
tions (Porte-Agel et al. 2001; Kleissl et al. 2003). The SGS 
Prandtl number (PrsGs) is more or less constant inside the 
boundary layer and gradually increases to <~ 1 in the inver- 
sion layer, as commonly assumed. Moreover, the values of 
Ptsgs increase in the surface layer. Earlier, we described 
very similar behavior in the case of turbulent Prandtl number 
(Pr t ). 



TABLE 3. z-less values of turbulence statistics. 
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FIG. 14. Correlation between u and w (r„ w ) (top), u and 6 (r u $) 
(middle), and w and 9 (r w $) (bottom). These statistics are averaged 
over the last one hour of simulation. 

5. Concluding Remarks and Future Perspectives 

One of the contributions of this research is the devel- 
opment and implementation of a new-generation subgrid- 
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FIG. 15. Vertical profiles of the SGS coefficients: C's, PrsGS, 
dynamically obtained by the locally-averaged scale-dependent dy- 
namic (LASDD) model. These profiles are averaged over the last 
one hour of simulation. 



scale model, termed as the locally-averaged scale-dependent 
dynamic (LASDD) model. This SGS model shares most 
of the desirable characteristics of the plane-averaged scale- 
dependent dynamic models, originally proposed by Porte- 
Agel et al. (2000) for the SGS stresses, and Porte-Agel (2004) 
for SGS fluxes. For example, it does not require any a priori 
specification of SGS coefficients since they are computed 
dynamically in a self-consistent manner. The dynamically 
estimated coefficients are found to strongly depend on fil- 
ter scale and atmospheric stability, in close agreement with 
a priori field studies (Porte-Agel et al. 2001; Kleissl et al. 
2003). However, in contrast to the original plane-averaged 
version, the LASDD model does not suffer from the in- 
sufficient SGS dissipation problem in simulations of stable 
boundary layers. 

The potential of our SGS model is made clear in coarse- 
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FIG. 16. Vertical profiles of the scale-dependent parameters: (3, 
Pe, dynamically obtained by the locally-averaged scale-dependent 
dynamic (LASDD) model. These profiles are averaged over the last 
one hour of simulation. 



resolution large-eddy simulations of moderately stable bound- 
ary layers. Overall, the agreements between our LES- 
generated turbulence statistics and observations, as well as 
some well-established empirical formulations and theoreti- 
cal predictions are remarkable. The results also show clear 
improvements over most of the traditional SGS models in 
the surface layer. In essence, we showed that tuning-free 
simulations of stable atmospheric boundary layers are fea- 
sible even with relatively coarse resolutions if one uses a 
robust and reliable SGS scheme. 

The next logical step would be to check the performance 
of this new-generation SGS scheme in simulating very sta- 
ble boundary layers. This would of course require exten- 
sive validation against existing profiles of various turbulence 
statistics measured during different field campaigns (e.g, 
Cooperative Atmosphere-Surface Exchange Study 1999 - 
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CASES99; Beaufort Sea Arctic Stratus Experiment - BASE; 
Surface Heat Budget of the Arctic Ocean - SHEBA). One 
of the characteristics of strongly stratified boundary layers 
is the existence of global intermittency (turbulent burstings 
in the midst of a laminar flow). In contrast to the traditional 
SGS models, the locally-averaged scale-dependent dynamic 
LES model has the correct behavior in laminar and transi- 
tional flows. This makes us believe that this model will be 
able to model the complex intermittency behavior of the very 
stable boundary layer flows. 

Contemporary SBL research has revealed that in very sta- 
ble regimes the relative importance of radiative cooling and 
heat exchange with the underlying soil becomes as signifi- 
cant as turbulent mixing (van de Wiel 2002). This means 
that radiation and soil physics should also be included in 
LES models before attempting very stable simulations. An- 
other interesting feature of this type of boundary layer flows 
is the presence of gravity waves. Conceptually, LES is 
capable of simulating gravity waves, provided the domain 
size is large enough. Unfortunately, the present computa- 
tional power dictates that in such cases one must have a rel- 
atively coarse resolution making thus the locally-averaged 
scale-dependent dynamic model more desirable than the 
laminarization-prone traditional SGS models. These issues 
will be addressed in our future research. 
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